/*
* INSANE - Interactive Structural Analysis Environment
*
* Copyright (C) 2003-2005
* Universidade Federal de Minas Gerais
* Escola de Engenharia
* Departamento de Engenharia de Estruturas
* 
* Author's email :    insane@dees.ufmg.br
* Author's website :  http://www.dees.ufmg.br/insane
* 
* This program is free software; you can redistribute it and/or
* modify it under the terms of the GNU General Public License
* as published by the Free Software Foundation; either version 2
* of the License, or any later version.
* 
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
* GNU General Public License for more details.
* 
* You should have received a copy of the GNU General Public License
* along with this program; if not, write to the Free Software
* Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
*/
package br.ufmg.dees.insane.model.femModel.problemDriver;

import java.io.Serializable;

import br.ufmg.dees.insane.model.femModel.element.Element;
import br.ufmg.dees.insane.model.femModel.element.ThinPlateElement;
import br.ufmg.dees.insane.util.IMatrix;
import br.ufmg.dees.insane.util.IVector;

public class VeubekeProblemDriver extends KirchhoffThinPlate implements Serializable {

	private static final long serialVersionUID = 1L;

	/* (non-Javadoc)
	 * @see br.ufmg.dees.insane.model.femModel.problemDriver.ProblemDriver#getC(br.ufmg.dees.insane.model.femModel.element.Element)
	 */
	public IMatrix getC(Element e) {
	    IMatrix cartesianNodalCoordsMatrix = ((ThinPlateElement)e).getCartesianNodalCoordsMatrix();
	    IVector c0 = this.getCoordsNodeCenter(cartesianNodalCoordsMatrix);
	    int numberOfDegreesOfFreedom = ((ThinPlateElement)e).getNumberOfDegreesOfFreedom();
		IMatrix finalStiffnessMatrix = new IMatrix(numberOfDegreesOfFreedom,numberOfDegreesOfFreedom);
	    finalStiffnessMatrix.setZero();
		IMatrix StiffnessMatrix = new IMatrix(numberOfDegreesOfFreedom,numberOfDegreesOfFreedom);
	    StiffnessMatrix.setZero();

	    double t = ((ThinPlateElement)e).getThickness();
	    double E = ((ThinPlateElement)e).getYoungsModulus();
	    double v = ((ThinPlateElement)e).getPoissonsRatio();

	    double D = E*t*t*t/(12*(1-(v*v)));

	    double a = Math.sqrt(Math.pow(cartesianNodalCoordsMatrix.getElement(0,0)-c0.getElement(0),2)+Math.pow(cartesianNodalCoordsMatrix.getElement(0,1)-c0.getElement(1),2));
	    double b = Math.sqrt(Math.pow(cartesianNodalCoordsMatrix.getElement(1,0)-c0.getElement(0),2)+Math.pow(cartesianNodalCoordsMatrix.getElement(1,1)-c0.getElement(1),2));
	    double c = Math.sqrt(Math.pow(cartesianNodalCoordsMatrix.getElement(2,0)-c0.getElement(0),2)+Math.pow(cartesianNodalCoordsMatrix.getElement(2,1)-c0.getElement(1),2));
	    double d = Math.sqrt(Math.pow(cartesianNodalCoordsMatrix.getElement(3,0)-c0.getElement(0),2)+Math.pow(cartesianNodalCoordsMatrix.getElement(3,1)-c0.getElement(1),2));
	    
	    IMatrix H = this.getH(cartesianNodalCoordsMatrix,v);
	    
	    IMatrix K1 = this.getK(this.getA(a,b),this.getB(a,b),this.getC(a,b),H);
	    IMatrix K2 = this.getK(this.getA(-c,b),this.getB(-c,b),this.getC(-c,b),H);
	    IMatrix K3 = this.getK(this.getA(a,-d),this.getB(a,-d),this.getC(a,-d),H);
	    IMatrix K4 = this.getK(this.getA(-c,-d),this.getB(-c,-d),this.getC(-c,-d),H);
	    
	    double cte1 = D/(6*a*a*a*b*b*b);
	    double cte2 = D/(6*c*c*c*b*b*b);
	    double cte3 = D/(6*a*a*a*d*d*d);
	    double cte4 = D/(6*c*c*c*d*d*d);
	    
	    K1.setScale(cte1);
	    K2.setScale(cte2);
	    K3.setScale(cte3);
	    K4.setScale(cte4);
	    
	    IMatrix J = new IMatrix(19,19);		J.setZero();
	    
	    J.add(this.getTransformationKi(K1,this.L1()));
	    J.add(this.getTransformationKi(K2,this.L2()));
	    J.add(this.getTransformationKi(K3,this.L3()));
	    J.add(this.getTransformationKi(K4,this.L4()));
	    
	    StiffnessMatrix.add(this.getFinalKi(J,getN(a,b,c,d)));
	    finalStiffnessMatrix.add(this.getFinalKi(StiffnessMatrix,this.getP(cartesianNodalCoordsMatrix,a,b,c,d)));
	    return(this.correctionStiffnessMatrix(finalStiffnessMatrix));
	}
	
	private IMatrix getA(double x, double y){
		IMatrix a = new IMatrix(3,10);
		a.setZero();
		
		a.setElement(0,1,(-x*y*y));   a.setElement(0,4,x*y*y);
		
		a.setElement(1,2,(-x*y*y));   a.setElement(1,5,x*y*y);
		
		a.setElement(2,1,x*x*x);      a.setElement(2,2,(-2*x*x*y)); a.setElement(2,3,(-6*x*x)); 
		a.setElement(2,4,x*x*x);      a.setElement(2,5,(-2*x*x*y)); a.setElement(2,6,6*x*x);
		a.setElement(2,8,(-2*x*x*y)); a.setElement(2,9,4*x*x*x);
		
		return a;
	}
	
	private IMatrix getB(double x, double y){
		IMatrix b = new IMatrix(3,10);
		b.setZero();
		
		b.setElement(0,1,(-2*x*y*y));  b.setElement(0,2,y*y*y);    b.setElement(0,4,(-x*y*y));
		b.setElement(0,5,(-y*y*y));    b.setElement(0,7,(-x*y*y)); b.setElement(0,9,4*x*y*y);
		
		b.setElement(1,1,(-x*x*y));    b.setElement(1,7,x*x*y);
		
		b.setElement(2,2,(-x*x*y));    b.setElement(2,8,x*x*y);
		
		return b;
	}
	
	private IMatrix getC(double x, double y){
		IMatrix c = new IMatrix(3,10);
		c.setZero();
		
		c.setElement(0,0,6*y*y);     c.setElement(0,1,x*y*y);    c.setElement(0,2,y*y*y);
		c.setElement(0,3,(-6*y*y));  c.setElement(0,4,2*x*y*y);  c.setElement(0,5,(-y*y*y));
		c.setElement(0,7,(-x*y*y));  c.setElement(0,9,4*x*y*y);
		
		c.setElement(1,0,6*x*y);     c.setElement(1,1,x*x*y);    c.setElement(1,2,x*y*y);
		c.setElement(1,3,(-6*x*y));  c.setElement(1,4,x*x*y);    c.setElement(1,5,(-x*y*y));
		c.setElement(1,9,4*x*x*y);
		
		c.setElement(2,0,6*x*x);     c.setElement(2,1,x*x*x);    c.setElement(2,2,x*x*y);
		c.setElement(2,3,(-6*x*x));  c.setElement(2,4,x*x*x);    c.setElement(2,5,(-2*x*x*y));
		c.setElement(2,8,x*x*y);     c.setElement(2,9,4*x*x*x);
		
		return c;
	}
	
	private IMatrix getH(IMatrix cN, double ni){
		IMatrix h = new IMatrix(3,3);
		h.setZero();
		double cos = this.getCosAlpha(cN);
		double sin = Math.sqrt(1-(cos*cos));
		double cte = 1/(sin*sin*sin);
		
		h.setElement(0,0,1);                    h.setElement(0,1,(-2*cos));                     h.setElement(0,2,(cos*cos+ni*sin*sin));
		h.setElement(1,0,(-2*cos));             h.setElement(1,1,(4*cos*cos+2*(1-ni)*sin*sin)); h.setElement(1,2,(-2*cos));
		h.setElement(2,0,(cos*cos+ni*sin*sin)); h.setElement(2,1,(-2*cos));                     h.setElement(2,2,1);
		
		h.setScale(cte);
		
		return h;
	}
	
	private double getCosAlpha(IMatrix cN){
		IVector c = this.getCoordsNodeCenter(cN);
		
		double cos = 0.0;
		double r1x = cN.getElement(0,0)-c.getElement(0);
		double r1y = cN.getElement(0,1)-c.getElement(1);
		double r2x = cN.getElement(1,0)-c.getElement(0);
		double r2y = cN.getElement(1,1)-c.getElement(1);
		
		cos = (r1x*r2x+r1y*r2y)/(Math.sqrt((r1x*r1x)+(r1y*r1y))*Math.sqrt((r2x*r2x)+(r2y*r2y)));
		
		return cos;
	}
	
	private IVector getCoordsNodeCenter(IMatrix cN){
		IVector c = new IVector(2);
		double x1 = cN.getElement(0,0);		double y1 = cN.getElement(0,1);
		double x2 = cN.getElement(1,0);		double y2 = cN.getElement(1,1);
		double x3 = cN.getElement(2,0);		double y3 = cN.getElement(2,1);
		double x4 = cN.getElement(3,0);		double y4 = cN.getElement(3,1);
		
		double x0 = (-1)*(x2*x3*y4 - x2*x1*y4 + x2*y3*x1 - x2*x3*y1 - x4*y3*x1 + x4*x3*y1 + x4*x1*y2 - x4*x3*y2)/(y1*x2 - y1*x4 - y3*x2 + y3*x4 - x1*y2 + x1*y4 + x3*y2 - x3*y4);
		double y0 = (-y1*x4*y2 + y1*y4*x2 + y3*x4*y2 - y3*y4*x2 + x3*y1*y2 - x3*y1*y4 - y3*x1*y2 + y3*x1*y4)/(y1*x2 - y1*x4 - y3*x2 + y3*x4 - x1*y2 + x1*y4 + x3*y2 - x3*y4);
		
		c.setElement(0,x0);		c.setElement(1,y0);
		
		return c;
	}
	
	private IMatrix getK(IMatrix A, IMatrix B, IMatrix C, IMatrix H){
		IMatrix K = new IMatrix(10,10);		K.setZero();
		
		IMatrix KA = new IMatrix(10,10);	KA.setZero();
		IMatrix KB = new IMatrix(10,10);	KB.setZero();
		IMatrix KC = new IMatrix(10,10);	KC.setZero();
		
		IMatrix auxKA = new IMatrix(10,3);	auxKA.setZero();
		IMatrix auxKB = new IMatrix(10,3);	auxKB.setZero();
		IMatrix auxKC = new IMatrix(10,3);	auxKC.setZero();
		
		auxKA.mulTransposeLeft(A,H);
		KA.mul(auxKA,A);
		
		auxKB.mulTransposeLeft(B,H);
		KB.mul(auxKB,B);
		
		auxKC.mulTransposeLeft(C,H);
		KC.mul(auxKC,C);
		
		K.add(KA);
		K.add(KB);
		K.add(KC);
		
		return K;
	}

	private IMatrix L1(){
		IMatrix l = new IMatrix(10,19);		l.setZero();
		
		l.setElement(0,0,1);
		l.setElement(1,1,1);
		l.setElement(2,2,1);
		l.setElement(3,3,1);
		l.setElement(4,4,1);
		l.setElement(5,5,1);
		l.setElement(6,6,1);
		l.setElement(7,7,1);
		l.setElement(8,8,1);
		l.setElement(9,15,1);
		
		return l;
	}

	private IMatrix L2(){
		IMatrix l = new IMatrix(10,19);		l.setZero();
		
		l.setElement(0,0,1);
		l.setElement(1,1,1);
		l.setElement(2,2,1);
		l.setElement(3,9,1);
		l.setElement(4,10,1);
		l.setElement(5,11,1);
		l.setElement(6,6,1);
		l.setElement(7,7,1);
		l.setElement(8,8,1);
		l.setElement(9,16,1);
		
		return l;
	}
	
	private IMatrix L3(){
		IMatrix l = new IMatrix(10,19);		l.setZero();
		
		l.setElement(0,0,1);
		l.setElement(1,1,1);
		l.setElement(2,2,1);
		l.setElement(3,3,1);
		l.setElement(4,4,1);
		l.setElement(5,5,1);
		l.setElement(6,12,1);
		l.setElement(7,13,1);
		l.setElement(8,14,1);
		l.setElement(9,18,1);
		
		return l;
	}

	private IMatrix L4(){
		IMatrix l = new IMatrix(10,19);		l.setZero();
		
		l.setElement(0,0,1);
		l.setElement(1,1,1);
		l.setElement(2,2,1);
		l.setElement(3,9,1);
		l.setElement(4,10,1);
		l.setElement(5,11,1);
		l.setElement(6,12,1);
		l.setElement(7,13,1);
		l.setElement(8,14,1);
		l.setElement(9,17,1);
		
		return l;
	}

	private IMatrix getTransformationKi(IMatrix Ki, IMatrix L){
		IMatrix K = new IMatrix(19,19);		K.setZero();
		IMatrix auxK = new IMatrix(19,10);	auxK.setZero();
		
		auxK.mulTransposeLeft(L,Ki);
		K.mul(auxK,L);
		
		return K;
	}

	private IMatrix getFinalKi(IMatrix A, IMatrix B){
		IMatrix K = new IMatrix(B.getNumCol(),B.getNumCol());		K.setZero();
		IMatrix auxK = new IMatrix(B.getNumCol(),B.getNumRow());	auxK.setZero();
		
		auxK.mulTransposeLeft(B,A);
		K.mul(auxK,B);
		
		return K;
	}
	
	private IMatrix getN(double a, double b, double c, double d){
		IMatrix N = new IMatrix(19,16);
		N.setZero();
		
		double n00 = c/(a+c);
		double n01 = -(a*c)/(6*(a+c));
		double n06 = a/(a+c);
		double n07 = (a*c)/(6*(a+c));
		double n012 = -(2*a*c*d)/(3*(a+c)*(b+d));
		double n013 = (2*a*c*d)/(3*(a+c)*(b+d));
		double n014 = (2*a*c*b)/(3*(a+c)*(b+d));
		double n015 = -(2*a*c*b)/(3*(a+c)*(b+d));
		
		double n10 = 3/(a+c);
		double n11 = -a/(2*(a+c));
		double n14 = d/(2*(b+d));
		double n16 = -3/(a+c);
		double n17 = -c/(2*(a+c));
		double n110 = b/(2*(b+d));
		double n112 = -(2*a*d)/((a+c)*(b+d));
		double n113 = -(2*c*d)/((a+c)*(b+d));
		double n114 = -(2*b*c)/((a+c)*(b+d));
		double n115 = -(2*a*b)/((a+c)*(b+d));
		
		double n22 = c/(a+c);
		double n28 = a/(a+c);
		double n212 = -(2*a*c)/((a+c)*(b+d));
		double n213 = (2*a*c)/((a+c)*(b+d));
		double n214 = -(2*a*c)/((a+c)*(b+d));
		double n215 = (2*a*c)/((a+c)*(b+d));

		N.setElement(0,0,n00);   N.setElement(0,1,n01);   N.setElement(0,6,n06);   N.setElement(0,7,n07);
		N.setElement(0,12,n012); N.setElement(0,13,n013); N.setElement(0,14,n014); N.setElement(0,15,n015);

		N.setElement(1,0,n10);   N.setElement(1,1,n11);   N.setElement(1,4,n14);   N.setElement(1,6,n16);
		N.setElement(1,7,n17);   N.setElement(1,10,n110); N.setElement(1,12,n112); N.setElement(1,13,n113);
		N.setElement(1,14,n114); N.setElement(1,15,n115);
		
		N.setElement(2,2,n22);   N.setElement(2,8,n28);   N.setElement(2,12,n212); N.setElement(2,13,n213);
		N.setElement(2,14,n214); N.setElement(2,15,n215);

		for (int i = 0; i < 16; i++) {
			N.setElement((3+i),i,1);
		}
		
		return N;
	}
	
	private IMatrix getP(IMatrix cN, double a, double b, double c, double d){
		IMatrix p = new IMatrix(16,16);
		IVector c0 = this.getCoordsNodeCenter(cN);
		
		double x0 = c0.getElement(0);		double y0 = c0.getElement(1);
		double x1 = cN.getElement(0,0);		double y1 = cN.getElement(0,1);
		double x2 = cN.getElement(1,0);		double y2 = cN.getElement(1,1);
		double x3 = cN.getElement(2,0);		double y3 = cN.getElement(2,1);
		double x4 = cN.getElement(3,0);		double y4 = cN.getElement(3,1);
		
		double x10 = x1-x0;			double y10 = y1-y0; 
		double x20 = x2-x0;			double y20 = y2-y0; 
		double x30 = x3-x0;			double y30 = y3-y0; 
		double x40 = x4-x0;			double y40 = y4-y0; 

		double x12 = x1-x2;			double y12 = y1-y2; 
		double x23 = x2-x3;			double y23 = y2-y3; 
		double x34 = x3-x4;			double y34 = y3-y4; 
		double x41 = x4-x1;			double y41 = y4-y1;
		
		double k11 = (x10*y20-y10*x20)/(Math.sqrt((x10*x10+y10*y10)*(x12*x12+y12*y12)));
		double k21 = (x20*y30-y20*x30)/(Math.sqrt((x20*x20+y20*y20)*(x23*x23+y23*y23)));
		double k31 = (x30*y40-y30*x40)/(Math.sqrt((x30*x30+y30*y30)*(x34*x34+y34*y34)));
		double k41 = (x40*y10-y40*x10)/(Math.sqrt((x40*x40+y40*y40)*(x41*x41+y41*y41)));
		
		double k12 = -(x10*x12+y10*y12)/(4*(x12*x12+y12*y12)*Math.sqrt(x10*x10+y10*y10));
		double k22 = -(x20*x23+y20*y23)/(4*(x23*x23+y23*y23)*Math.sqrt(x20*x20+y20*y20));
		double k32 = -(x30*x34+y30*y34)/(4*(x34*x34+y34*y34)*Math.sqrt(x30*x30+y30*y30));
		double k42 = -(x40*x41+y40*y41)/(4*(x41*x41+y41*y41)*Math.sqrt(x40*x40+y40*y40));
		
		double Xx = (x1-x3)/(a+c);
		double Yx = (y1-y3)/(a+c);
		double Xy = (x2-x4)/(b+d);
		double Yy = (y2-y4)/(b+d);
		
		p.setElement(1,1,Yx);  p.setElement(1,2,Xx);
		p.setElement(2,1,Yy);  p.setElement(2,2,Xy);

		p.setElement(4,4,Yx);  p.setElement(4,5,Xx);
		p.setElement(5,4,Yy);  p.setElement(5,5,Xy);

		p.setElement(7,7,Yx);  p.setElement(7,8,Xx);
		p.setElement(8,7,Yy);  p.setElement(8,8,Xy);

		p.setElement(10,10,Yx);  p.setElement(10,11,Xx);
		p.setElement(11,10,Yy);  p.setElement(11,11,Xy);
		
		p.setElement(12,0,-6*k12);    p.setElement(12,1,y12*k12);   p.setElement(12,2,x12*k12);  p.setElement(12,3,6*k12);
		p.setElement(12,4,y12*k12);   p.setElement(12,5,x12*k12);   p.setElement(12,12,k11);
		
		p.setElement(13,3,-6*k22);    p.setElement(13,4,y23*k22);   p.setElement(13,5,x23*k22);  p.setElement(13,6,6*k22);
		p.setElement(13,7,y23*k22);   p.setElement(13,8,x23*k22);   p.setElement(13,13,k21);
		
		p.setElement(14,6,-6*k32);    p.setElement(14,7,y34*k32);   p.setElement(14,8,x34*k32);  p.setElement(14,9,6*k32);
		p.setElement(14,10,y34*k32);  p.setElement(14,11,x34*k32);  p.setElement(14,14,k31);
		
		p.setElement(15,0,6*k42);     p.setElement(15,1,y41*k42);   p.setElement(15,2,x41*k42);  p.setElement(15,9,-6*k42);
		p.setElement(15,10,y41*k42);  p.setElement(15,11,x41*k42);  p.setElement(15,15,k41);
		
		return p;
	}
	
	private IMatrix correctionStiffnessMatrix(IMatrix k){
		IMatrix t = new IMatrix(16,16);
		t.setZero();
		t.setElement(0,0,1);
		t.setElement(1,2,1);
		t.setElement(2,1,-1);
		t.setElement(3,3,1);
		t.setElement(4,5,1);
		t.setElement(5,4,-1);
		t.setElement(6,6,1);
		t.setElement(7,8,1);
		t.setElement(8,7,-1);
		t.setElement(9,9,1);
		t.setElement(10,11,1);
		t.setElement(11,10,-1);
		t.setElement(12,12,1);
		t.setElement(13,13,1);
		t.setElement(14,14,1);
		t.setElement(15,15,1);
		
		return this.getFinalKi(k,t);
	}

	public String getUnknownDualVariableLabel() {
		// TODO Auto-generated method stub
		return null;
	}
}
